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ABSTRACT 


Computer simulation studies of two frequency domain 
adaptive beamforming algorithms for planar arrays are 
presented. The algorithms are modified complex LMS adaptive 
algorithms that can process an arbitrary number of 
harmonics. The algorithms provide estimates of the spherical 
coordinates (1.e., range, depression angle, and bearing 
angle) of multiple broadband targets in both the near-field 
and far-field. Computer simulation results comparing the 
average estimation error for range, depression angle, and 
bearing angle as a function of the input SNR, range (near- 
field and far-field), and harmonic number, are presented. 
The "full angular coverage" capability of the algorithms was 


also tested. 
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I. INTRODUCTION 


Frequency domain adaptive beamforming is a_e signal 
processing technique used to cophase the output electrical 
Signals at each element in an array of sensors. The adaptive 
algorithm recursively adjusts the complex weights at each 
element in the array minimizing the mean-square-error between 
a reference signal and its estimate. The resulting steady- 
state phase weights represent the uncorrupted phase variation 
across the face of the array. The phase values at the 
elements in the array contain the location information of the 
source. Herein lies the usefulness of such an algorithm. By 
cophasing the output signals at each of the elements in a 
receive array uSing an adaptive beamformer, noise corruption 
that occurs in the medium and/or at the receiver can be 
effectively cancelled. The resulting steady-state phase 
values can be manipulated to determine the position of the 
source relative to the array. One obvious application for a 
beamformer of this nature is in determining the position of 
a sonar target in an ocean mediun. 

Target localization in bearing and depression angles has 
been successfully achieved using an adaptive beamformer of 
this design [Refs. 1, 2]. However, one localization parameter 
that is of significant interest and yet unsolved for by the 


previously mentioned signal processing algorithm [Refs. 1, 2] 


is "target range." The objective of this thesis is to develop 
a localization algorithm that can process spherical wavefront 
information in a noise environment in order to determine the 
position of a target in range, bearing, and elevation. The 
frequency domain adaptive beamformer developed in [Refs. 1, 
2] will be applied to this problem. 

In our analysis, we modeled the target as a broadband 
sound source. As we convert a received output electrical 
Signal to the frequency domain, the complex frequency spectrum 
of the target contains many frequency components. The 
frequency domain adaptive beamforming algorithm processes each 
spectral line of the target's signature independently. Asa 
result, the target's position is estimated for each frequency 
component. Given a multi-target situation, the algorithm is 
capable of localizing each target provided that at least one 
unique spectral line can be associated with each sound source 
(CargeT). 

In Chapter II, we will develop the theoretical groundwork 
used in solving the localization problem. First, we consider 
the acoustical properties of wave propagation in an isospeed 
ocean medium. In Section II.A, we begin with a general form 
of the inhomogeneous wave equation and develop the output 
electrical signals at each of the elements ina planar array. 
Based on a given signal-to-noise ratio (SNR), we corrupt the 
output electrical signals with noise. This time-domain signal 


is converted to the frequency domain for processing (Sec. 


eigen The signal is filtered using a frequency domain, 
modified least-mean-squares (LMS) adaptive algorithm. Since 
we are analyzing spherical waves, this routine has been 
modified to include analysis that assumes non-separable 
complex weights. Other phase distortions introduced by the 
signal processing routine are cancelled using techniques 
discussed in the remainder of Section II.B. The last section 
of Chapter II presents an algorithm that was developed to 
process the complex weights in order to determine the position 
of the target in spherical coordinates. 

In Chapter III, several computer simulation studies are 
discussed. The test cases were selected in order to 
demonstrate the significant properties of the localization 
algorithm. In each case, the run was first conducted in a 


noise-free environment to validate the propagation and 


Mecalization models. Gwee sethe wbaselime — results ) were 
generated, the test cases were repeated in a noise 
environment. 


Section III.A tests the full angular coverage capability 
of the algorithm. Full angular coverage implies that the 
target can be localized regardless of its relative position 
to the receive array. A single target is placed at a 
broadside and then endfire position relative to the array. 
The range of the target is then varied to determine the 
performance of the algorithm in the near and far-field 


regions. One additional test case placing the target at an 


arbitrary geometry relative to the array is also examined in 
this section. In Section III.B, the multi-harmonic capability 
of the algorithm is exercised. The test case reviews the 
performance of the algorithm given a single target with 
several harmonics, covering a wide frequency spectrum. In the 
final section of our results chapter, Section III.C, we 
examine the multi-target performance of the algorithm. In 
this simulation study, three targets with different locations 
are tested. The final conclusions and recommended areas for 


further research are presented in Chapter IV. 


II, THEORETICAL DEVELOPMENT 


The purpose of this chapter is to present the theoretical 
groundwork used in solving the localization problem. The 
approach can be divided into three major areas (see Figure 
yee) . In generating the output electrical signals at the 
elements in the receive array, we first consider the 
acoustical nature of the problem. Beginning with a general 
form of the wave equation, we derive a mathematical model for 
spherical wave propagation in a homogeneous medium. The wave 
propagation model is used to generate both the input acoustic 
Signals and the output electrical signals at the transducers 


in the array. Based on a given signal-to-noise ratio (SNR), 
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Figure 2.1 The Block Diagram Representation of the Analysis 
Steps Used in Solving the Localization Problem. 


the output electrical signals are corrupted by noise. Using 
a frequency-domain adaptive beamformer, the output signals 
from the elements in the receive array are co-phased. The 
steady-state phase weights determined by the beamformer are 
then processed by a localization algorithm that provides an 


estimate of the target's position in spherical coordinates. 


A. SIGNAL GENERATION 

The starting point of the localization problem is the 
generation of the acoustic signal incident on the planar 
array. The signal must represent wave propagation in an 
isospeed ocean medium. Specifically, our signal must exhibit 
spherical wavefront curvature that is a function of the range 
between the source and the array. This mathematical signal 
model is developed by solving the wave equation. The 
corresponding output electrical signals are then discretized 
by taking a preset number of time samples over the data record 
length at each element. The final step of signal generation 
1s to add noise to the electrical signals. The output from 
this section is a time-domain, noise-corrupted signal at each 
transducer element in the receive array (see Figure 2.2). 

1. Spherical Wave Propagation Model 

The propagation of acoustic energy in the ocean can 


be described by the following inhomogeneous wave equation: 
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Figure 2.2. The Block Diagram Representation of the Major 
Steps Performed in Signal Generation. 


where o(t,r) is the velocity potential at time t and position 
r, Xy(t,r}) represents the source distribution, and c(r) is the 
speed of sound in the ocean. In our problem, we are 
interested in wave propagation in a homogeneous medium. 
Therefore, c(r) 1s set equal to a constant sound speed c. 
When the speed of sound is constant, the solution to the wave 


Equatiten (22 ))) issq1ven by- 


me ( (So) 


2a 
|r -ro| ( 


] 
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Vo 


where | r-Fo | is shown in Figure 2.3 [R@£.93-p. Zoo 

Next, we model the source distribution generated by 
the target as an omnidirectional point source with arbitrary 
time dependence located at r = r, (see Figure 2.3). That is, 


let 


xm(tr) = g(t) 6(r-ro) (2.3) 


where g(t) is an arbitrary function of time. By substituting 
Equation (2.3) into Equation (2.2) and performing jen 


integration, the velocity potential g(t,r), reduces to 


Gr, 
| oes (2.4) 
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Figure 2.3 Spherical Coordinates I= | rol, 89, and, Wo defining 
the Position of the Target. 


where 


R=|r-Po |. (2.5) 


As a result, we can express the acoustic signal Yu(t, vr) 


incident upon the array [Ref. 4] as 


(es) 
ym(t.r)= ao . (22> 
At this point, we can apply the physical conditions 
of our array to this general expression of the input acoustic 
Slgnal Equation (2.6). First, we shall assume that the output 
electrical signal y(t,r) is directly proportional to the inpme 


acoustic signal, that is, 
Vil ean Vin ey (2.90) 


where the constant of proportionality has been dropped. Next, 
we W1ll sample this continuous electrical signal at a sampling 
rate f,., providing a discrete form of the time-domain signal 
in preparation for signal processing. Consider a planar array 
lying in the XY plane. Given the int2relement spacings of the 
array dad, and d,, the element position vector r (see Figure 


2.3) can be expressed as 


r=md,x + ndyy, 5 (0). (2.65 


10 


The location of the source as measured from the center of the 


Gienay x, 1s 


ro = XoX + Yoy + Zo Z. (2.9) 


Therefore, the magnitude of the range to the target, R, can 


be expressed as 


(ee 10) 


yp 
ie 


R=| r-Fo |=| (md, - x9)? +(ndy- yo) + 20. 


By taking L time samples over the record length, the 


discrete representation of the signal is 


qm | ) 
= C27 
y( IT., md,, ndy } = 5 
where "1" corresponds to a particular time instant and 
T =! 
ae oe (2.12) 


is the sampling period. 
The arbitrary function of time g(t) is represented by 


the following finite Fourier series: 


K 
e(t)=ap+2), a,cos ( 2nqfot +0, ), (ols) 
oo 


sual 


where: a, is the DC component of the source signal, 

a, 1s the magnitude of the Fourier series coefficient 
of the source signal at harmonic gq, 

©. is the phase of the Fourier series coefficient of 
the source signal at harmonic gq, 

f. is the fundamental frequency of the source signal, 


and K is the total number of harmonics. 


By setting time t in Equation (2.13) equal to the retarded 

time [ir,-8) and substituting this equation into Equation 
Cc 

(2.6), we obtain the final form of the output electrical 


Signal at the receive array: 


K 
y( 9) Gera 9 eae don anh agcos| 2mafo(tr,-® a | - (2.0% 
q= 1 


2. Noise Addition 

To validate the signal processing and localizatzon 
algorithms, as well as our wave propagation model used in 
generating the signals, we first perform the analysis ina 
noise-free environment. However, once the baseline results 
have been obtained, we must test the performance of the 
overall algorithm under realistic environmental conditions. 
Therefore, we corrupt the output signal at each element in 
the array with noise. Additive, zero mean, white, Gaussian 


noise 1S added to the output signals from each element in the 


Ny 


planar array. We select various signal-to-noise ratio (SNR) 
values to simulate a particular noise environment. The 
resulting signals represent discrete time, noise corrupted, 


output signals at each element in the planar array. 


B. FREQUENCY DOMAIN SIGNAL PROCESSING 

The output signal at each of the transducer elements is 
converted to the frequency domain by using a Discrete Fourier 
Transform (DFT) routine. In this section, we will begin our 
discussion by deriving the form of the radiated signal after 
taking the DFT. Next, we will examine all phase components 
that surface in our frequency-domain analysis. Common terms 
associated with wavefront curvature will then be presented. 

The frequency-domain Signal processing routine can be 
Givided into three major functions (see Figure 2.4). The 
modified least-mean-sgquares (LMS) adaptive algorithm will be 
the first topic of discussion. By cophasing the output 
Signals from the elements in the array to yield a minimum 
least-squares error, the noise corruption of the signal is 
reduced. Another source of phase distortion occurs in the 
conversion to the frequency domain. The phase values 
calculated by taking the DFT are limited to a closed interval 
between ([-7,7]}. This phase "wrapping" effect must be 
rectified prior to further processing by the localization 
algorithm. This effect is corrected by the phase unwrapping 


algorithm. The "unwrapped", steady-state phase weights are 
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Figure 2.4 Block Diagram Representation of the Major 
Steps Performed in Signal Processing. 
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then smoothed along the X and Y axes using non-linear least- 
squares estimation. The output from this signal processing 
section is an estimate of the uncorrupted phase at each 


element in the receive array. 


is Phase Components Encountered in Frequency Domain 
Analysis 


From our analysis in Section II.A.1, we developed a 
discrete, time-domain expression of the output electrical 
Signal at each of the elements in the planar array. AS in 
most Signal processing applications, it is often more 
convenient to operate on the signal in the frequency domain. 
Therefore, our first manipulation is to take the DFT of the 


eeemat. The sampled output spectrum, Y.(q,m,n), 18 given by 
~ j2nql/L 
me(4,m,n) = Cy, » ys tiatele nd, }e cans qm. K. C25) 


where eee lee alle 2 (2a Ga) 


eee 15 the complex weight at each element (m,n). 


SeeecacucLing Equation (2.11) into Equation (2.15) yields 


aaa Cee eK (2.17) 


7 
Y,(q.m,n) = +> » g(1r,-8 


ales 


The Fourier series coefficients introduced to represent the 
arbitrary function of time g(t) in Equation (29973) (Gane 


expressed by the following complex coefficients, c,: 


c=a.e % "Q=]l)2ee (2.18) 


Whence an is the magnitude and 0, is the phase of harmonic q. 


Since the Fourier series coefficients c, represent g(t), then 


we know from the properties of Fourier series that g(t + 


e) 


will have the Fourier Series Coettrerenes Cc ¢ + ies 


q 
Therefore, Equation (2.17) can be written in the following 


Simplified form (Ref. 5:pp. 162-166) 


I 3 
Y.(q,m,n) = pcre ea loRie = 1, ee (2 Sip 
where R/c is substituted for €,. The range R is given by 
Equation (2.5) as the magnitude of the vector ,»which is 


the range from a particular element in the array to the target 
(see Figure 2.3). Using vector algebra, this magnitude can 


be rewritten as 


} 


|r-ro|=[(r-ro}*(r-ro] (2.20) 


By expanding the dot product within the radical, k jeanmiee 


rewritten as 
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1/2 


2 a 2 
Jr-r0 [=| ro -2ro(a,,er |+r @22521) 





where a, = Ugx + Voy + Woz 62722325 


is the unit vector in the direction of r,, and is defined in 


mBemmsSeof direction» cosines u,, Vv and® was tet lowSea(see 


Oo? 


mreure 2.3): 


Ug=sin8g cos Wo (22) 
Vo=sinO9 sin Wo (220) 
and Wo = cos Oo. (225 ) 


By taking the dot product of Equation (2.8) with Equation 
(2.22), and substituting the result into Equation (2.21) we 


obtain: 
,) , 1/2 
|r -r9|= Tg - 2ro(ugmd, + vond,) + (md,} + (nd) ; (2.26) 
From this range expression and by referring back to Equation 


(2.19), we can define the phase term due to the propagation 


of acoustic energy from the source as 





6.=- : To - 2rq(ugmd, + vond,) + (md,)° + (nd, 7 (2.27) 


iby, 


This term, 6 is the phase due to source radiation and 


re? 


contains the physical localization information of the target. 


By solving for r,, u,, and v,, we can determine the spheraeam 


o! 
coordinates of the source. 

At this point in our discussion it is proper to 
introduce all of the other phase components that exist at a 
given transducer element in the receive array. Beginning at 
the source, the first phase term that we have discussed is 
due to the physical propagation of acoustic energy, 6. 
Additionally, another phase term exists that originates from 
the target. The source signal also contains an initial phase 
term at each harmonic. Specifically, we defined this in 
Equation (2.12) as the phase of the Fourier coefficient, O.- 
Since this phase is not known a priori from the localization 
Side of the problem, we must solve for this unknown in order 
to obtain the radiation phase term, 9.. Our objective in this 


section is to develop an approximation of the uncorrupted 


Signal phase and to pass this output 


0.= 9,49, (2.28) 


to the localization algorithm (see Section II.C). 
Three other phase terms exist that distort the signal 


phase mentioned above. We can represent the noise corruption 
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SEBCNeeInNpUL aGe@ustic Signal by awphase term, G+.@eThis phase 
term varies at each element and for each harmonic analyzed. 
The objective of the LMS adaptive algorithm is to cancel this 
noise phase term at each element. Another source of phase 
distortion occurs in the conversion to the frequency domain. 
The phase values calculated by taking the DFT are limited to 
a closed interval between [-7,7]. When a phase value exceeds 
the interval of [-,7]}, the phase is "wrapped" within the 
limits of the interval by adding/subtracting 27. This phase 


"wrapping" effect can be accounted for by the term 0 The 


wrap ° 
function of the "phase unwrapping algorithm" is to unwrap the 
phases at the elements along the X and Y axes. Lastly, a 
constant phase term is added to the phase at each element in 
the array as a result of the cophasing process of the LMS 
adaptive algorithm. This phase term, 9,, is not cancelled by 
the signal processing routine of this section. The phase term 
© is passed along with the estimate of ©,;, to the 
localization algorithm. This algorithm groups the 6, term 
with the ®, component, and cancels them both yielding the 0, 


term. In summary, the total phase component at a particular 


element can be described by the following equation: 


ets Orato art Ont’ wrap cars (2.29) 


IL, 


2. Common Terms Associated with Wavefront Curvature 

A wavefront is defined as a surface of constant phase. 
Figure 2.5 depicts two sound sources (targets) located at 
broadside relative to a planar array. The target shown in 
Figure 2.5.a is located close to the receive array. Looking 
at the incident wavefront on the surface of the array, we 
observe significant wavefront curvature. When a wavefront, 
as measured at the receive array, exhibits significant 
curvature, the target is considered to be in the "near-field 
(NF)". Figure 2.5.b shows the same target at a more distant 
position from the array. The amount of wavefront curvature 
seen at the array face 1s noticeably less. If a wavefront has 
slight curvature or resembles a plane wave, the target is 
considered to be in the "far-field (FF)". One observation 
that can be made from this figure is that the range 
information of a target's position can be found in the 
curvature of the wavefront. The range to the physical 
boundary between the near-field and far-field regions is a 
function of the array size and the wavelength of the incident 


acoustic field (see Figure 2.6) and is given by [Ref. 5:p. 33] 


DE T R array (253109 
rX 
where p is the range to the NF/FF boundary, R is the 


array 


effective radius of the array, and X is the wavelength of 
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Figure 2.5 Representation of the Range Effects 
on Wavefront Curvature. 
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Figure 2.6 Representation of Near-Field/ Far-Field Boundary. 
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mmemineident £f£1eld- By considering a rectangular array of M 


x N elements, Equation (2.30) becomes 


p=B (bp) +(4y4) 


ee be 


(2221) 

where 
eg Silo Lai. (82) 
ey SIN (238) 
d,= y= ain, (2.34) 

and 

denin= Spa (2735) 
Given a target with K,, harmonics, f,,, 1s equal to 


mwew total number of harmonics, K,.,, times the fundamental 
frequency, f., that is 


0 


max 10: ( 22536 ) 


Since the broadband signal of the target can encompass a wide 
range of frequencies, one unique value for p that delineates 
the boundary between the near and far fields does not exist. 
Rather, one can only specify a unique NF/FF boundary for a 
given harmonic of interest due to the frequency dependence of 
the range, p. However, a target can be classified as a near- 
field target at all harmonics provided that it is located 
within a minimum p value defined in terms of the maximum 


wavelength as follows: 


Pao: 


er UAlAl ean 


Conversely, if a target's range is greater than a maximum 9) 
value, it is classified as a far-field contact for all 


harmonics, where 
Pou 7| (4) AL (bA)} (2 en 


This NF/FF boundary is a helpful tool in classifying the 
relative target position to the array since it takes into 
account array size and target frequency components. One 
misconception that has resulted from establishing a boundary 
of this nature is the argument that only targets located in 
the near-field can be localized in range. Although it is true 
that a near-field target has significant wavefront curvature, 
our research demonstrates that an adequate range estimate can 
be determined given a wave with a slight amount of wavefront 
curvature (i.e., for targets located well into the far-field). 

At the beginning of this section, we defined a 
wavefront as a surface of constant phase. Using the phase 
information at the elements of the planar array, another way 
in which wavefront curvature can be defined is by plotting 
the phase as a function of element position across the X and 
Y axes. For a near-field target, this phase distribution will 


exhibit significant curvature. The phase distribution across 
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Beemer the X Gr Y axis of the array for a far-field target 
will be more linear. In the results section, we will plot the 
ideal phase distribution across these two array axes for each 
of the cases that we analyze. 
3. Least-Mean-Squares (LMS) Adaptive Algorithm 

As discussed in the previous section, the phase passed 
to the frequency-domain modified LMS adaptive algorithm for 
processing 1S comprised of many components. The measured 


phase, ©, can be described by 


m? 


0 n= 8,+ 95+ Ont Ourap- Zo) 


The objective of the adaptive algorithm is to cancel out the 
feece term, ©... The algorithm performs this function by 
adaptively updating the complex weights at each element in 
the array while attempting to minimize the mean-square error. 

The LMS adaptive algorithm is a well known signal 
processing aid [Ref. 6]. In recent research work [Ref. 1], 
a frequency-domain modified LMS adaptive algorithm was applied 
to a planar array to process plane-wave Signals. 
Specifically, complex weights at each element in the array 
were optimized adaptively with respect to a reference signal. 
The ideal value of the phase of each complex weight is equal 
to the negative value of the uncorrupted wrapped signal phase. 


The algorithm used for this plane wave case assumed that the 


complex weights were separable. This 1s a valid assumption 


ELD) 


for plane-wave analysis or cases when the wavefront curvature 
is minimal (i.e., far-field problems). However, for those 
cases with significant wavefront curvature and a resulting 
interdependence between x and y coordinates, this assumption 
is not valid. 

The following equations summarize the key steps that 
comprise the adaptive algorithm assuming separable complex 
weights [Ref. 2}. For the remainder of the thesis, this 
algorithm will be called LMS,,.,, where the word “flat” ainsi 


negligible curvature. 
“ ] 
Estimate: zila)=iypqd dilan)dc(4,m) Y{q,m,n } (2.40) 
n m 


where z;(q) is the estimate of the reference signal at harmonic 


h 


q at the i‘ iteration and where the complex weights c, are 


assumed to be separable, that is, 


Cmnl 4) = q,m)d{q,n} (2.41) 
Error Siemals: e.(q)=z(q)-z;(q) (2.42) 


Complex Weight Update: 


Ci+1(q,.m)=c;(q,m)+ 2u;(q)ei(q) > 4)(q.n} Y(q.m,n} (2.43a) 


n 


(2.43b) 
dis 1(q.n)=d;(q.n)+2u;(q)e;(q}} ¥ c;(q,m) Y(q,m,n} 
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L(q)= 


2 
where EMa0.0)] 24 3c) 


is the step size parameter. The adaptive routine computes 
mamations (2.40) through Equations (2.43c) N times, minimizing 
Eemerror signal e,. After each iteration, the magnitude of 
the complex phase weights c,,,(q,m) and d,,,(q,n) are normalized 
ioeorder to maintain unit magnitude. By normalizing the 
magnitude of the phase weights after each iteration, the phase 
component of the phase weights is the term that is recursively 
optimized to minimize least-squares error. With this 
modification, the final output from the algorithm is a set of 
steady-state, phase weights that represent the negative values 
of the uncorrupted wrapped phase of the output electrical 
Signal at each element in the array. 

The adaptive algorithm summarized above (derived in 
[Ref. 2]) was modified to process near-field waves (1.ée., 
waves with significant wavefront curvature) by expressing the 
complex weights in a non-separable form. This algorithm will 
Memeecterred to as IMS |. where the word "curve" implies 
Significant curvature. The complex weights c, are set equal 


to w(m,n) in these equations: 
“ l 
Estimate: 25110 = era yl CC Y{q.m,n } (2.44) 
m n 


Earer Signal: e(q)=z(q)-z(q] 45) 
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Complex Weight Update: 
Wi+1(4.m,n)=w;(q,m,n)+ 2u;(q)e)(q) Y (q.m,n} (2. Aen 


where uw.(q) is defined in Equation(2.43c). As in the LMS;,., 
algorithm, the magnitude component of the complex phase 
weights is normalized to unity after each iteration. 


In summary, adaptive algorithm LMS assumes that the 


flat 
complex weights are separable and will work best for cases 
with minimal wavefront curvature. LMS wurve assumes Nici 
separable weights and will give better performance for cases 
with significant wavefront curvature. In the results section, 
both algorithms have been used for each scenario presented to 
demonstrate these properties of the two algorithms. 
4. Phase Unwrapping Algorithm 

The output from the two LMS adaptive algorithms are 
processed sets of steady-state, "wrapped" phase values of the 
output electrical signals at each element in the array. As 
mentioned earlier, phase distortion occurs in the conversion 
to the frequency domain. The phases calculated by the DFT are 
limited to a closed interval between [-7,7]. As a phase value 
exceeds this interval, the phase is "wrapped" within the 


limits of [-7,7] by adding/subtracting 27. This phase 


"wrapping" effect is undesirable and must be rectified prior 
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to further processing by the localization algorithn. The 
phase unwrapping algorithm performs this task. 

Figure 2.7 is an example of the complexities of phase 
wrapping. Figure 2.7.a depicts the ideal, "unwrapped" phase 
Variation across the X axis. As the phase value exceeds +180° 
(see Figure 2.7.b), a noticeable jump in phase is evident as 
the phase is wrapped within the limits. This jump is detected 
by the phase unwrapping algorithm and is properly corrected. 

The phase unwrapping algorithm compares the phase 
difference between two adjacent elements. If the difference 
exceeds +m radians, the unwrapping algorithm is activated and 
a phase correction is made. The routine then moves to the 
next adjacent element and performs the same logic check and 
correction as necessary. Prior to discussing the details of 
the logic steps performed in the unwrap algorithm, we will 
first examine the basis for choosing +7 as the criteria for 
phase wrap detection. The maximum phase difference resulting 
from wave propagation (@,) that could exist between adjacent 
elements 1s egual to 7. If a phase difference exists that is 
greater than 7, then it must be due to phase wrapping. 

The two extreme cases that yield the maximum phase 
slope across the X axis are shown in Figure 2.8. The first 
example 1s a contact located in the extreme near-field at a 
broadside orientation to the array. An additional constraint 
is that the maximum slope only occurs at the highest harmonic. 


The second example 1s a target located at endfire relative to 
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Figure 2.7 The Complexities of Phase Wrapping. 
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Figure 2.8 Cases Where the Phase Difference Between Elements 
Approaches the Limiting Value of 180 Degrees. 


the array and positioned in the far-field. The maximum slope 
for this case also occurs at the highest harmonic only. 
Starting with our expression for the radiated phase Ouy(eee 
Equation (2.27)), we will apply the conditions of these two 
extreme cases and calculate the maximum slope. 

First, let us rewrite Equation (2.27) for use in this 


Section. 


2naf 1/2 


2 
oe =| To - 2ro(ugmd, + vond,) + (md,)° + (nd,)” (2.495 





For each case we will only perform the analysis across the X 
axis, since similar results can be derived assuming the Y 
axis. Therefore, we set the index in the Y direction, n= 0O. 
In the limiting near-field case, the target is located at 
broadside relative to the array. As a result, the direction 
cosine term u, is equal to zero. For a target positichecmaas 
the very near field, the square of the range | tei 
considered negligible as compared to (md,)°*. Therefore, 


Equation (2.47) can be simplified as follows: 


4 on Ke agai Nal. 
ers € 


— 


(2.48) 


where K.,, is the highest harmonic. The derivative of OG yjyageg 


respect to element number m is equal to 


be _ 2m | pacpyas W((Ss 
max ; 


dm % (2.49) 
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Paewsing | baviatronm (2.54) through Equation (2.36), the 


interelement spacing term d, can be expressed as 
lon (2.50) 


Substituting Equation (2.50) into Equation (2.49) yields 


dO, 


eam = + (ono 1a) 
ae 


For the second case, the target 1s located at endfire. 
Therefore, the maximum value that u, can attain 1s +1. If we 


eereeu, — +1, then Equation (2.47) simplifies to 


1/2 
Jaf) ote 
_hK max'0 ro? ~ 2ro(md, )+ (md, (2752) 


—r 
— 


: Cc 


mae term inside the radical of Equation (2.52) can be 


rewritten as 
2 eee 2 2 
ko. = To| x }+(md,] = (To - md, | : (2953) 
Therefore, the radiated phase can be expressed as 


9 a4 TK mato | | ne (2.54) 


Taking the derivative with respect to m yields 
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Ohi) GAG 
Kies = pao (2.55) 
dm max Cc 

Finally, substituting Equati¢6mo 250) nem ae into Equaties 


(2.55) yields the maximum slope 


e =e ie (2.56) 
dm / max 

The phase unwrapping algorithm uses this value of 7 
as the criterion for detecting a phase wrap between elements. 
The logic flow of the phase unwrapping algorithm is pictured 
in Figure 2.9. The phase unwrapping algorithm begins at the 
center element of the array. If the unwrapping is to be 
performed in the (+) X direction, then the step size is set 
equal to unity (for negative X direction, s = -1). Two 
iterations through the comparator section of the routine are 
required. This double pass approach is necessary due to the 
effects of wavefront curvature. That is to say, since the 
Slope of the phase variation can be both positive and negative 
across the axis, the direction of unwrapping becomes more 
complex. In the first pass, the phase values between adjacent 
elements are compared. When a wrap is detected a positive 
correction is applied to the phase value. This positive 
correction 1s equal to +2ka, where the constant k 12S9igme 
number of times the phase value at the adjacent element has 


been wrapped. On the second pass through this routine, a 
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Start at Center 
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If unwrap along 
+ X axis 






Two unwrap 
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Do Count = 1,2 





Figure 2.9 Logic Flow of the Unwrapping Algorithm. 
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backup check of the phase unwrapping iS conducted. If the 
phase value has been incorrectly unwrapped, the algorithm 
applies a negative correction. In this way, even in cases of 
Significant wavefront curvature, the values are correctly 
unwrapped between the elements. 

This process is repeated along the (+) X axis out to 
the final element. The routine then returns to the center 
element and iterates along the (-) X axis. Upon unwrapping 
the X axis, the algorithm operates on the Y axis in a similar 
Manner. Since the center element is a common element for both 
axes, it 1s essential to start the unwrapping here for both 
cases (versus beginning at the end of either axis). The 
output from this section is a set of unwrapped, steady-state, 
phase weights at the elements along the X and Y axes. The 
final step in the signal processing routine is to smooth the 
phase values along the X and Y axeS uSing non-linear least- 
Squares estimation. 

5. Non-Linear Least-Squares Estimation 

The two major sources of phase distortion, that is, 
6, and 6,,.,, have been eliminated using the techniques 
mentioned in the last two sections. One further step is 
necessary prior to passing the unwrapped, steady-state phase 
information to the localization algorithm. Although the 
modified LMS adaptive algorithms cancel the majority of the 
noise corruption in the received signals, some jitter still 


exists in the steady-state phases. The localization algorithm 
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uses the phase information at three elements along the X and 
Y axes to estimate the range to the target. Specifically, the 
phase values at the end elements and the center element are 
used. Our research has shown that the slight jitter that may 
exist in the phase values at these elements can have a radical 
impact on the resulting range estimate. Therefore, the last 
step in our Signal processing routine is to find a curve that 
represents the best fit given the discrete phase values along 
the two axes. The technique used to perform this smoothing 
1S non-linear least-squares estimation. 

Given n+l equally spaced data points assigned to 
values of x, where x = 0, 1, .. . =, n, a non-linear curve 
that best fits this data can be determined using orthogonal 
polynomials [Ref. 7]. The equation of this curve is given by 


(Ref. 7] 


P(x} = ag Pio(x} + ay Pa (x) Baal Tela ee Bg 2 57) 


where m is the degree of the polynomial, 1 lS ara cr a) 


Seyi x) are the orthogonal polynomials, and a,, a,,...) a, are 
the coefficients of the orthogonal polynomials. Orthogonal 
polynomials have the following property: 

(2.58) 


y Px) Pnx(x) = 9, j#k. 


x=0 
The general formula used in computing the orthogonal 


polynomials is [Ref. 7] 


ou 


m+ | (x)' 
’ got) J Bees n 
, : 2.53 
(n)' ( ) 


Pam (x)= Len? I”; 





where ( 2yn6iGe 


Heo 


Using this formula, the first four polynomials are given as 





follows: 
P io(x} = 1, (2 688) 
Pailt}=1-2-, (2.620 
eae og Xie- (2.63) 
Sa EO n(n- 1)’ 
and 
x-] x - 1}(x - 2] 
P 3(x) = 1 - 12— +305 i 205 (2.64) 


By using these four polynomials and by calculating the values 
for their respective coefficients, we can construct a cubic 
polynomial to fit the given phase data while minimizing the 
least-squares error. The equation for calculating the 


coefficients is given by [Ref. 7} 


S> fix) Pasta 
=0 


xX 
a. = ———_————._. 1= 0, 1,...,m, (2.655 


S) P2 x) 


c=0 
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where the values of f(x) are the phase values given for each 
element along a particular axis. 

Using orthogonal polynomials, one can derive the 
equation for a curve of degree m with minimum least-squares 
emror, where m < n (n+l is the number of data points). For 
the results that we will present in Chapter III, we performed 
a best fit approximation using a cubic polynomial (m = 3). 
A polynomial of this degree was adequate to perform the 
smoothing required for our phase data. 

One additional value that we calculated and monitored 


in our computer runs was the actual least-squares error. This 


error can be calculated from the following equation [Ref. 7]: 


m 


n n 

B= > fix - Yad, Pri]. (2.66) 
xe) 1=0 x=0 

The value of this error gives an indication of how well the 

non-linear curve fits the data. 

This concludes all the steps used in the frequency- 
domain signal processing section of our research. The output 
from this section is a set of smoothed phase values for the 
elements along the X and Y axes where the effects of noise 
have been minimized and the frequency-domain signal has been 
smoothed, and unwrapped. The target localization algorithm 
Will take these values and determine an estimation of the 


spherical location of the target. 
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C. TARGET LOCALIZATION IN SPHERICAL COORDINATES 

Beginning with the noise-corrupted output electrical 
signals at the elements in the receive array, our goal has 
been to take this time-domain information and localize a 
target or a group of targets in spherical coordinates. One 
intermediary step was to convert the signals into the 
frequency domain and filter out as much of the distortion due 
to noise as possible. With this accomplished, we are now in 
a position to take these processed signals and work backwards 
in order localize the targets. 

The localization algorithm applies the physical properties 
of wave propagation developed in Section II.A, and utilizes 
the symmetry of the receive array to obtain a position 
estimate. The first step in this routine is to isolate the 
terms that contribute to wavefront curvature, and then to 
estimate the target range. Having solved for this parameter, 
the wavefront curvature is cancelled, leaving a linear phase 
relationship. By manipulating the slope information from this 
conditioned phase data, the angular position of the target is 
determined (see Figure 2.10). 

1. Three Point Range Algorithm 

Our first step in the localization problem is to 
manipulate the unwrapped, smoothed, steady-state phase 
information received from the signal processing routine to 


estimate target range. We denote this processed phase as 0.- 
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Figure 2.10 Block Diagram Representation of the Major Steps 
Performed in the Localization Algorithm. 
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As discussed in section II.B.1, 6, is comprised of three phase 


components, that is, 


0,=0,+6,+0, (2.67) 


where © is the phase due to wave propagation, Qo, is the 
phase of the Fourier series coefficient c,, and 0, is the phase 
term that represents a constant phase shift that is generated 
at each element due to the LMS algorithm. In our range 


algorithn, @, and ©, are grouped together and considered as one 


constant phase term Og that is, 
Oy = 6,4 6.. (2.68) 


Recall from Equation (2.27) that, the radiated phase is given 


by 
1/2 


6.(q,m,n)= - 9 - 2rp(ugmd, + vond,) + (md,}’ + (nd,}” , (2.69) 





27Gfo 
Cc 


The range algorithm uses the phase information at 
three elements per X and Y axis to estimate range. For this 
analysis, let us consider the X axis (set n=0). The phase of 
the two end elements on each side of the X axis are combined 
to form one equation in terms of range r, and 6,,. Since the 
end elements are symmetrical about the center element, this 
equation can be simplified. The second equation is formed by 


the value of the phase at the center element. This equation 
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is also a function of our two unknowns, Yr, and 6,,. These two 
equations are manipulated to solve for the two unknowns. In 
this fashion, we obtain a value for the range based on three 
elements along the X axis. The same approach is implemented 
on the three corresponding elements on the Y axis. AS a 
result, we obtain two estimates of the range. 

The processed phase at the center element is obtained 
by combining Equation (2.67) through Equation (2.69) and 


setting m=0 and n=0: 





2taQf 9 
9 ,(q,0,0)=0,.- z 


To- (27.570) 
This is the first equation needed to solve for our two 
unknowns. The equation from the symmetrical end elements 
requires a more detailed derivation. First, from Equation 


(2.67) and Equation (2.68), we can write that 
9 (q,m,n)— 6 ,:= 6,(q,m,n}. (Poll) 


Assuming M elements along the X axis, the element number 


total 


m of these two symmetrical elements can be expressed as +m', 


where m' is defined as 


m' = (Miotai- 1)/2. (2372) 
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Let us derive an equation at element -m' first. To eliminate 
the radical of Equation (2.69), we square both sides of 
Equation (2.71) after substituting Equation (2.69) into 
Equation (2.71). The element number -m' is then substituted 


in place of m yielding 


a? 
2. | 2 4n q fol 2 wae 
9, (a-m‘.0) — 26 {a-m\0) 6 |. + 0g 20g ug{-m'}d, +(-m df. (2.7m 


Cc 


In the same manner, the equation at +m' is given by 


2 
An q’ fo 


2 
C 


2 2 
0 ,{a.m'.0) — 28, (a.m,0) 6 1+ 6, := |r. 2rgugm'd, +(m' af) : (2 745) 


Equation (2.73) 1s added to Equation (2.74) to cancel out the 


term containing the direction cosine information to yield 


2 2 2 8x q fof2_, ,, 2) (2a 
0 ,(a.-m'.0) + 8 ,{a.m‘,0) — 26418 pla-m.0) + 8 (am) + 20,:= —— +g + (m'd,} | S ) 


This is the second equation needed to solve for our two 


unknowns. Next, we square both sides of Equation (2.70) to 
£OrM 
Pee dy? 
2 2 4n q fo 2 
,0,0 aa Vy . —_——— ae ° 
8, (9.0.0) — 2614.0.) 8 ,-+ 8, oe (2.76) 


This equation is then substituted into Equation (2.75) to 


solve for the 6. estimate, 


44 


2 2 2 8n a7 f2 7 
O; 
205 — On Op at —z d,} 


a 2[20,,-(2,..+ pa] 


which is one of the two unknowns. Having solved for the first 


(2n7 7) 


@) 
i 


? 


unknown, we substitute the estimated value for 84: IMtco 
Equation (2.70) to solve for r,. Upon algebraically solving 


for r, in terms of 6,, from Equation (2.70), we obtain 


Fo=[0q~ 8pla.0.0)][ 5], (2.78) 

In summary, we cancelled the term that contains the 
Girection cosine variable u, in Equation (2.69) by using 
symmetrical elements along the X axis (tm'). We combined the 
resulting equations derived from these two elements with the 
equation representing the phase at the center element. From 
these three elements, we formed two equations to solve for two 
unknowns. One of these unknowns is the range estimate to the 
target, Ty. The routine is then repeated for the 
corresponding three elements along the Y axis. 

2. Estimation of Bearing and Depression Angles 

Knowing the range estimate % and the value for 0. 
greatly simplifies the task of solving for the angular 
estimates. We manipulate the terms containing these two 
quantities within the phase at each element and effectively 
cancel out the wavefront curvature along the X and Y axes. 


We perform a linear least-squares estimation to calculate a 
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line that best fits this corrected phase data. The slope of 
this line is equal to the direction cosine u, given data along 
the X axis; v, for Y axis data. From the direction Cocum@e. 
estimates, we can calculate the bearing and depression angles. 

For this angular analysis, our goal is to formulate 
a new phase term that is a function of a single direction 
cosine only. Specifically, we solve for the term u,m by 
setting n=O in Equation (2.69) and define this value as our 
new phase term. The first step is to subtract the estimated 
value of O,: found in the previous section from the phase value 
at each element, to solve for © (see Equation 2.71). We 
substitute this value for 6, into Equation (2.69) and solve 
for Um as follows 


2 
c ~ 2 2 
anak, On ign (mee 


-2r9d, 


By setting m=O and performing a similar algebraic 


manipulation, v,n can be solved for as 


0, (q,m,n} = vgn = - , m=O. 


Using these two equations, we recalculate the new phase at 
each element along the X and Y axes. These conditioned phases 
are linear along the X and Y axes. From the non-linear 
estimation discussion of Section II.B.5, we utilize the first 


two orthogonal polynomials P,, and P,, (see Equation (2.61) and 
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Equation (2.62)) to perform a linear least-squares estimation 
on this data. The slope of the resulting line estimates is 
Peliciato Uisandev, Gespecelvely. 

Been ene Semical cCulatedmVvalles mor UU. .and v,, swe can 
calculate the estimates of the bearing and depression angles 


as follows: 


1/2 
O5= sin” |( ud+v3] | (2.81) 


and 


“~ -1{ Yo 
Wo = tan Pal (2582) 


Pethe this information, our localization of the target is 


complete (see Figure 2.3). 
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III. SIMULATION RESULTS 


Three major categories were chosen in order to validate 
the significant capabilities of the localization algorithm. 
In Section III.A, the results validating the full angular 
coverage capability are presented. Full angular coverage 
implies that the target can be localized regardless of its 
relative position to the receive array. The multi-harmonic 
capability of the algorithm is exercised in Section III.B. 
The two test cases presented review the performance of the 
algorithm given a single target with several harmonics, 
covering a wide frequency spectrum. In the final section of 
this chapter, Section III.C, we examine the multiple broadband 
target performance of the algorithm. Three targets with two 
unigue harmonics each and at different locations are tested. 

Prior to presenting the results from these test cases, a 
review of the design limitations of our computer simulation 
1s necessary. The computer code is written in the Fortran 
language and iS organized in modular sections (1i1.e., 
subroutines reflect major analysis blocks in Figures 2.4 and 
2.9). Since the program was written as a research tool, the 
emphasis was placed on making the code detailed and easy to 
follow. One adverse consequence of this decision is longer 


processing times. Given the restrictions on computer 
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simulation run times, certain self-imposed design limitations 
were required. 

The first design limitation is the number of time samples 
taken when converting the signal to the frequency domain. [In 
all of the test cases conducted, 65 time samples at each 
element were taken. In the context of our research, we found 
that our localization estimates improved when we took more 
time-domain samples over the data record. This principle is 
also formally documented in ([Ref. 1]. The other major 
constraint was the number of elements that comprised the 
planar array. In each of the following test cases, we used 
an 11 x 11 element array. Two beneficial results occur as we 


increase the number of elements in the receive array. First, 


the dimensions of the array increase. With a larger array, 
the wavefront curvature becomes more pronounced. Secondly, 
the performance of the modified, frequency-domain, LMS 


adaptive algorithm, as well as the curve fitting accomplished 
by the non-linear least-squares estimation routine are 
enhanced as the data from more elements is processed. 

In each of the test cases, the simulation run was first 
conducted in a noise-free environment to validate the 
propagation and localization models. Once the baseline 
results were generated, the test caseS were repeated in a 
noise environment. The two noise environments’ chosen 
represented SNR values of 0.0 @B and -3.0 @B. An overview of 


the algorithm's performance is presented at the beginning of 
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the respective test case. The average estimation error of the 
target's location at a given harmonic are presented for each 
test case. 

The modified, frequency-domain LMS adaptive algorithms 
conduct a specified number of iterations as the complex 
weights are recursively updated (see Section II.B.3). In all 
of the noise cases that were tested, the number of iterations 
used in the LMS adaptive algorithms was equal to 100. The 
estimation errors that are presented for each test case are 
average values. For a given SNR, the average estimation 
errors for the range, bearing and depression angles were 
obtained by running the computer simulation 50 times. 

Lastly, since the range estimate 1s very sensitive to 
wavefront curvature, a plot of the ideal phase variation 
across the X or Y axis is given for each test case. Aside 
from the actual range position (1i.e., NF/FF), other factors 
that affect wavefront curvature include frequency (li.e., 
harmonic number) and relative geometry. Rather than examine 
the contribution of each of these individual effects, the 
performance of the range estimate will be considered in light 
of the plot of wavefront curvature for a given target 


(location) and harmonic number. 
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A. CATEGORY I: VALIDATION OF THE "FULL ANGULAR COVERAGE" 
CAPABILITY 


A single target was placed at three positions relative to 
the center of the receive array. The first geometry analyzed 
is with the target located at a position broadside relative 
to the receive array. The range of the target is then varied 
to examine the range limitations of the algorithn. In the 
second test case, an arbitrary position between broadside and 
endfire geometries is analyzed. For these two cases, a single 
frequency component of 1000 Hz was used in analyzing the 
target. In the final case, the target is placed at a position 
endfire to the plane of the array. Three harmonics with a 
fundamental frequency of 1000 Hz were used for the target in 


mies case. 


1. Case I.A: Target Located Broadside Relative to the 
Planar Array 


In this case, the target is located at a depression 
mere O| Of O10" (1.e., broadside relative to the array). In 
Table 1, the no-noise data is presented. Three range cases 
are summarized in this table. At a range of 5.89 meters, the 
target is located well within the near-field (0.1*p). A 
second range located at a position in the middle of the NF 
region, 29.4 meters (0.5*p) was tested. The third range 
examined represents a position in the extreme far-field, 58.9 
km (1000%*p) and was chosen to show how well the algorithm 


works in this ideal no-noise case, given minimal wavefront 


oak 


TABLE 1 
CASE 1LA.1- TARGET AT BROADSIDE, SINGLE HARMONIC, NO NOISE. 


Actual Location Average Estimation Error Range % Diff 


ee e e e 
Case (m) (deg) (deg) (m) (m) (deg) (deg) . 
Oulr5 
Flat ' 
Curve 


Flat 
Curve 


Flat 
Curve 





D2 


eurevacure. The plot of the wavefront curvature for these 
three range values is given in Figure 3.1. 

The results from the two modified LMS adaptive 
algorithms are presented in Table 1. The algorithm that 


assumes separable weights, LMS,,,,, is identified as "Flat" in 


the table. The algorithm assuming non-sSeparable phase 
memonts, EMS 4, 1s denoted by "Curve". Recall that LMS,,,, 
is ideally suited for minimal wavefront curvature. This 


property is supported by the results in that, as target range 
increases, the range estimation error decreases. In fact, 
when the target is located in the extreme NF region, the range 
estimate using LMS,,,, 1s 7% larger than the actual range 
value. iio epersern=a. well aerall three of these range 
cases. 

Note the wavefront curvature in Figure 3.1 for the FF 
case (1000*p). For this no-noise case, both algorithms 
accurately estimated this extreme range value. For all the 
results presented, both algorithms correctly estimated the 
mepGessiOn angle 6,. In this case, the bearing angle is 
arbitrary. Throughout the remaining examples, we will only 
discuss the highlights of the results. Since both algorithms 
generally give good estimates of the angular information, most 
of our discussion will be focused on range estimates. 

Table 2 summarizes the broadside results for the 0.0 
dB SNR case. The LMS..,. algorithm does well in the near 


field. However, as the range position is increased to a value 


a 
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TABLE 2 
CASE J.A.2- TARGET AT BROADSIDE, SINGLE HARMONIC, 0.0 DB. 


Actual Location Average Estimation Error Range % Diff 
e c 

Bo to] ER, “Ry 
(deg) _ (deg) ‘i 


-Q.1 
-0.3 


























-174.6 
-166.2 










-Q.1 
-0.3 


-168.9 
-172.8 










-0.1 
-0.3 


-167.2 
-164.8 








-0.1  -163.0 
-0.3 -180.7 









-0.1 
-0.2 


-185.0 
-192.9 





ae 


into the FF, the range estimates become unreliable (i.e., 
range estimation error is large; X and Y estimates differ) due 
to minimal wavefront curvature. LMS,,,, estimates range out to 
a position of 3*p adequately. Figure 3.2 shows the difference 
in wavefront curvature for 3*p to 5*p, the region where the 
range estimates of the LMS,,,, algorithm become erroneous. In 
the remainder of our test cases, the range estimates of TENaam 
will be shown as superior to those of LMS... as we enter the 


FF region. However, LMS... will be seen to perform better in 
the extreme NF. This property reflects the assumptions made 
concerning the complex phase weights (1.e., separable versus 
non-separable). Even as range estimates became 
unsatisfactory, both algorithms continued to estimate the 
angular information to within a few tenths of a degree. 

Table 3 summarizes the broadside results for -3.0 @B 
SNR. AS expected, in this noisier environment, the estimates 
of both algorithms are degraded. Once again, since range is 
the most sensitive of the estimated spherical coordinates, the 
corresponding estimates are more severely degraded. Note that 
the 6, estimate is still accurate even as the range esta. 
becomes more unreliable. 

2. Case I.B: Target with Arbitrary Position 

For this test case, the target is positioned at an 
arbitrary position given a bearing angle ¥, of 140° and Wa 
depression angle ©, of 35°. Figure 3.3 depicts the wave tees 


curvature for the range values tested in this case. At this 
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TABLE 3 
CASE IL.A.3- TARGET AT BROADSIDE, SINGLE HARMONIC, -3.0 DB. 


Actual Location Average Estimation Error Range % Diff 
e e e 
ums{| Ro 9% % R, Ry 80 “YW, CR, eR, 
Case (m) (deg) (deg) (m) (m) (deg) (deg) - 


Flat zQ 5 Sioa 
0.5 -167.2 


-0.1 -168.5 
-0.4 -172.0 


-O.1 -174.7 
-0.4 -164.0 


-0.1 -162.8 
-0.4 -187.5 


-0.1 -184.9 
-0.3 -193.2 


-0.1 -176.5 
-0.4 -174.8 
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geometry, the wavefront curvature is considerably less than 
the broadside case. The resulting range estimates in a noise 
environment bear witness to this fact. 

Table 4 documents the results of the no-noise case. 
The results indicate minimal error even at a range of 1000*p. 
Tables 5 and 6 summarize the performance of the two algorithms 
in a noise environment. The same trends that were discussed 
in the last section are evident. Once again, LMS,,,, gives the 
most reliable maximum range (about 3*p). The angular 
estimates are excellent in both noise environments. 

3. Case I.C: Target with Endfire Geometry 

This test case is significant in that it exposes two 
limitations in the target localization problem. One 
limitation occurs in the phase unwrapping routine performed 
in the signal processing section. The second deficiency 
observed in this test case can be traced to the range 
estimation equation in the localization algorithm 
Specifically, at this extreme endfire position, the estimates 
of the localization algorithm become unreliable. 

The first discrepancy occurs when analyzing the 
highest harmonic of a target located at an endfire position 
relative to the X or Y axis. From our discussion in Section 
II.B.4, we identified this extreme case as the basis for 
choosing m radians as the setpoint for the phase unwrapping 
algorithm. In our research, we observed that with noise 


corruption, or with the slight distortion that resulted from 
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TABLE 4 


CASE I.B.1- TARGET AT ARBITRARY POSITION, SINGLE HARMONIC, 
NO NOISE. 


Actual Location Average Estimation Error Range % Diff 


a | 
Rp 8M | Rn Ry “Op WH ] OR, ER 
Case (m) (deg) (deg) (m) (m) (deg) (deg) 

Flat 
Curve 


Flat 0.03 0.02 
Curve 0.0 0.0 


Flat 
Curve 
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TABLE 5 


CASE I.B.2- TARGET AT ARBITRARY POSITION, 
SINGLE HARMONIC, 0.0 DB. 


Actual Location Average Estimation Error Range % Diff 


e eC ce € 
LMS Ro 9, Y Ro. Ry oT Wo CR “Ra 
Case (m) (deg) (deg) (m) (m) (deg) (deg) ‘i 
Flat 
Curve 


Flat 
Curve 
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TABLE 6 


CASE I.B.3- TARGET AT ARBITRARY POSITION, 
SINGLE HARMONIC, -3.0 DB. 


Actual Location Average Estimation Error Range % Diff 


ec e C 
rms| Ro %  % Re Ry 80 Mo} PR, fR, 
Case (m) (deg) (deg) (m) (m) (deg) (deg) x 


161 1972 





Os 


the signal processing section of our analysis (in the absence 
of noise), the phase difference between two elements might 
deviate slightly above the maximum slope of the phase 
variation (m radians). The phase unwrapping algorithm 
classifies this phase difference aS a wrap, even though one 
aid Not ececuE. The localization algorithm cannot properly 
process the phase data with this erroneous unwrapping. 
Conversely, the phase unwrapping algorithm might not detect 
an element where the phase had been wrapped. This unwrap 
anomaly only occurs at the highest harmonic when the target 
is located at endfire on an axis of the receive array. 

The second limitation observed from this endfire case 
1S more general in that it applies to all harmonics. 
Additionally, the discrepancy occurs when the target is 
located at a position endfire to the array, irrespective of 
the on-axis condition. Referring to Equation 2.77, which is 
one of the two equations essential for estimating the range, 
it can be seen that as the denominator approaches zero, the 
eguation becomes invalid. Therefore, one limit of this 
equation is when twice the phase value at the center element, 
0,(q,0,0), 1s equal to the sum of the phase values at the end 
elements, 0©,(q,-m',0) + 06,(q,m',0). This condition iss 
approached when the target is located at a position endfire 
to the receive array. 

Since the second limitation is more general than the 


unwrap anomaly, we chose to examine a target at an endfire 
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position off the array axis. Specifically, the bearing angle 
teeenosen 1s equal to 237°. Additionally, we examined a 
target with three harmonics with a fundamental frequency of 
1000 Hz. From Table 7, the no-noise case, we observe a 
discrepancy in the spherical coordinate estimates for LMS,,., 


in the extreme NF (0.1*p Tables 8 and 9 document the 


min? : 
performance of the algorithm in the presence of noise. Aside 


from the 0.1*p,.. case for LMS,,,,, both algorithms estimate the 


min 
angular information satisfactorily. However, the range 
estimates of both algorithms using the phase values in the Y 
direction are erroneous. The phase variation curves shown in 
Figures 3.4 and 3.5 were generated for the highest frequency 
component (3000 Hz). By comparing the curvature in Figure 3.4 
with that shown in Figure 3.5, we see that the phase variation 
along the Y axis is more linear. With near linear phase 
variation at endfire, the denominator of Equation 2.77 
approaches zero. 

In summary, this range estimation limitation only 
occurs at endfire geometries when a near-linear phase 
variation exists. The error introduced to the 84, estimate 
(see Equation 2.75) is cancelled by the resulting error in 
range estimation (see Equation 2.76) such that the angular 
estimates are still valid (see Equations 2.79 and 2.80). This 
discrepancy combined with the unwrap deficiency outlined above 
are the only two limitations observed as we analyzed the full 


angular coverage capability of the algorithm. 
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TABLE 7 


CASE I.C.1- TARGET AT ENDFIRE POSITION, 
THREE HARMONICS, NO NOISE. 


Average Estimation Error Range % Diff 


e e 
ES (m) (deg) _(deg) m wy de ) fea Ro. Ni 
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TABLE 8 


CASE I.C.2- TARGET AT ENDFIRE POSITION, 
THREE HARMONICS, 0.0 DB. 


— Estimation Error Range % Diff 


“6 
OE ae) oe - 


124 -0.1 
1.41 0.105 
0.861 0.06 
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TABLE 9 


CASE I.C.3- TARGET AT ENDFIRE POSITION, 
THREE HARMONICS, -3.0 DB. 


Actual Location Range % Dif 
e e€ e 
(m) (deg) (deg) m ; Wy dep deg Ros Roy 
Q.1* min 


-3.52 -26.0 
-0.15 -22.8 
-34.1. *-O:92 
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B. CATEGORY II: VALIDATION OF THE MULTI-HARMONIC CAPABILITY 

The wavefront curvature, as measured at the surface of a 
receive array, is a function of a target's range (see Section 
TI.B.2). In addition, the amount of wavefront curvature is 
also a function of the wavelength and thereby, frequency of 
the harmonic being processed. Given a target at a fixed 
position in space, for a given harmonic, the range value of 
the target can be represented by a constant times the range 
to the NF/FF boundary (const * p). As we analyze a target at 
a different harmonic, the value for the range to the NF/FF 
boundary (fp) changes (see Equation 2.31). As we discuss the 
results of this section, we will occasionally refer to the 
range to the NF/FF boundary for a given harmonic as p,, where 
q is the harmonic of interest. To describe this harmonic 
dependency of p, and the resultant relative position of the 
target to this p, consider the following example. Suppose a 
target has 10 harmonics to be analyzed in the frequency 
domain. Assume that at the highest harmonic, the range of the 
target can be specified as 0.1 * p,,. The range to the same 
target 1S equal to 1.0 * p, at the lowest harmonic. At the 
highest harmonic, the relative range to the target can be 
considered in the extreme near-field region. While, for the 
lowest harmonic, the relative range of the target is equal to 
p,, the range to the NF/FF boundary. From this example, we 
see that a scenario could exist where a target may be 


considered as a NF target at one harmonic and a FF target at 


Val 


another harmonic, even though the target remains at a fixed 
pOSltion ingspacer 

In this case study, we analyze the performance of the two 
algorithms given a target that has many frequency components. 
In the first case, the target is located at a relative range 
of 0.56 * ps (NF) at the highestiiamic ie At the lowest 
harmonic, the relative position of the target is equal to 5 
* pf, (FP). The second exanple 1s chosen to represent a target 
with a relative range in the extreme NF at the highest 
component (0.1 * p,), to a relative range closer to the NF/FF 
boundary at the lowest component (0.7 * 9,). 

1. Case II.A: Target With Nine Harmonics 

In this test case, a fundamental frequency was chosen 

at 400 Hz. With nine harmonics, the band of the frequency 
spectrum varies from 400 to 3600 Hz. The results summarizing 


the performance of the two localization algorithms are 


presented in Tables 10 through 12. In the no-noise case 
(Table 10), the IMS... data reflects the dependence Ota 
relative target position on frequency. At the lowest 


harmonic, where the target is at five times the relative NF/FF 
boundary (5 * p,), the range estimation error is largest. As 
harmonic number increases, the relative position of the target 
with respect to the p value at a given harmonic decreases. 
At the highest harmonic, the relative position of the target 
is in the NF and the range estimation error is the lowest. 


This effect 1s not evideneein ene IMS,;,,. results. 
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TABLE 10 
CASE II.A.1- SINGLE TARGET, MULTI-HARMONIC (9), NO NOISE. 


Actual Location Average Estimation Error Range % Diff 
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TABLE I1 
CASE I1.A.2- SINGLE TARGET, MULTI-HARMONIC (9), 0.0 DB. 


Actual Location Average Estimation Error Range % Diff 
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-2.46 
-1.18 

-0.350 
-0.027 
-0.016 
-0.143 
-0.034 


7.02 
-2.73 
-1.28 
-0.813 
-0.010 
-0.094 
-0.225 
-0.245 
-0.289 
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-0.559 
-0.311 
0.010 
-Q.096 
-0.017 
0.032 
0.007 
-0.003 
0.017 


-4.72 
-1.82 
-2.17 
0.062 
-0.164 
-0.093 
-0.079 
-0.218 
-0.364 


3.00 
2.48 
-0.855 
-0.843 
0.524 
0.268 
-0.094 
0.829 
-0.104 


SiS) 
25-0 
6.11 
2.08 
2.10 
-0.056 
-0.165 
0.587 
0.605 





TABLE 12 
CASE II.A.3- SINGLE TARGET, MULTI-HARMONIC (9), -3.0 DB. 
Average Estimation Error Range % Diff 
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Figure 3.6 describes the phase variation across the 
X axis for three harmonics. At the highest harmonic, the 
wavefront curvature is most significant; whereas, for the 
lowest component, the curvature is minimal. The test case 
results for the scenarios with noise corruption are presented 
in Tables 11 and 12. In general, LMS,,,, performs the best 
over the wide range of frequencies. For the lowest three 
harmonics, neither algorithm provides an adequate range 
estimate. 

2. Case II.B: Target With Seven Harmonics 

In this case, we examine a target with a range value 
equal to 0.7 * p,. The fundamental frequency is 400 Hz and 
the band of the frequency spectrum covers 400 to 2800 Hz. 
Since the lowest harmonic describes the maximum relative range 
of the target to the NF/FF boundary, the relative ranges at 
every harmonic will be within the NF. This case was chosen 
to demonstrate the performance of both LMS algorithms at each 
harmonic. As was previously discussed, LMS,,,, should perform 
best in cases with minimal wavefront curvature. The converse 


1s true for LMS A different bearing and depression angle 


curve ° 
were chosen for this case. 

The trends discussed in the last section are evident 
in the results presented in Tables 13 through 15 for this 
case. In the noise cases (Tables 14 and 15), LMS,,., gives a 


better range estimate for the first two harmonics than the 


LMS estimates. LMS performs best for the higher 


curve curve 
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TABLE 13 
CASE II.B.1- SINGLE TARGET, MULTI-HARMONIC (7), NO NOISE. 


Actual Location — Estimation Error Range % Diff 


eC € eC 
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TABLE 14 
CASE II.B.2- SINGLE TARGET, MULTI-HARMONIC (7), 0.0 DB. 


Actual Location Average Estimation Error Range % Diff 
pa Ro 9, Yo | . Rox : Ry ap 0 “W, CR S R. 
) | (m) (deg) (deg) m m) deg deg & a 
7* 


lat | 0.7*pmin 
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TABLE 15 
CASE II.B.3- SINGLE TARGET, MULTI-HARMONIC (7), -3.0 DB. 


Actual Location Average Estimation Error Range % Diff 


ae a 
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harmonics. The phase variation for harmonics three through 
seven contain significant wavefront curvature as seen in 
Pergure 3.7. Harmonics one and two have minimal wavefront 
curvature. This case demonstrates the regions where one LMS 
algorithm is superior to the other. 
C. CATEGORY IIT: VALIDATION OF THE MULTIPLE TARGET 

CAPABILITY 

In this final section, the performance of the two 
algorithms is evaluated in a multiple target, multiple 
harmonic case. Three targets are oriented at random 
geometries in both the NF and FF regions. Each target 
radiates two unique spectral lines resulting in six total 
harmonics. A fundamental frequency of 100 Hz was chosen. 
Therefore, the band of the frequency spectrum is 100 to 700 
Iz. Many factors affect the amount of wavefront curvature 
for a given harmonic. This test case shows the complexities 
involved in analyzing a target when several of the wavefront 
curvature limiting effects are combined. Two major factors 
are evident as we analyze the multiple targets in this 
scenario. First, as a target's depression angle increases 
from the broadside position, the amount of wavefront curvature 
becomes less significant. Secondly, as we analyze data at 
lower harmonics, the wavefront curvature is also less. These 
two factors were examined individually in Sections III.A.3 and 


Section III.B. 
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marge: bois located in the near field (0.3 * p,..) at 
@ = 42° and ¥ = 153° and radiates harmonics 1 and 3. Target 


2 is closer to the NF/FF boundary (1.8 * p,..) at @, = 71° and 


min 
~ = 218° (harmonics 2 and 5). The final target is in the FF 
mo * pf) at 6, = 2° and ¥ = 47° (harmonics 4 and 6). In 
order to prevent the FF target from being masked by the noise 
levels of the other targets, the initial magnitude of the 
Fourier series coefficients were weighted so that the signal 
strength of each target at the receive array would be equal. 
We will summarize the results in Tables 16 through 18 by 
considering the targets one at a time. 

From Table 16, in the absence of noise, the LMS,,,,, routine 
localizes Target 1 with minimal estimation errors. However, 
Iisa noise environment (see Tables 17 and 18), IMS,,., yields 
a more accurate estimate of the spherical coordinates at both 
harmonics. Given the phase variation shown in Figure 3.8, it 
is clear that at both harmonics, significant wavefront 
Piature 1S present. Under these conditions, LMS... should 
perform better than LMS,,,,, but it does not. 

This performance discrepancy can be attributed in part to 
the design of the array. Since the interelement spacing is 


optimally designed for the highest frequency of interest, the 


performance of the localization algorithm at lower frequencies 


is degraded. The results recorded in Table 16 bear witness 
to this fact. The range estimation error is maximum at the 
lowest harmonic. A similar trend was also observed in the 
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TABLE 16 
CASE HI- MULTI-TARGET (3), MULTI-HARMONIC (6), NO NOISE. 


Average Estimation Error Range % Diff 
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TABLE 17 
CASE I0- MULTI-TARGET (3), MULTI-HARMONIC (6), 0.0 DB. 


Average Estimation Error Range % Diff 


| (m) __(deg) _(deg) m m deg) (deg) Ox by 


stale, «= -0275 «= 1.99 E29 
72.5  -18.8 -0.596 -0.10 
-0.542 0.051 = 1.58 1.46 

626 -53.1  -0.044 -0.270 
2.96 -2.52 -0.301 -0.047 
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TABLE 18 
CASE Ill- MULTI-TARGET (3), MULTI-HARMONIC (6), -3.0 DB. 


Actual Location Average Estimation Error Range % Diff 


0 0 0 y 
a (m)__(4eg) _ (ew Roy 


-0.488 1.85 be 
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Q:023 Gisz 1.49 
145.7. -0.070 -0.241 

-4.91 -0.280 -0.038 


12700 = =0.050 = 0.055 





86 


‘L# JOBIe | *[]] IseD AoJ 


€ aIn3iq 


b s1uowe PY JO UOIPUN v SB SIXe-X JY} SSOIDY UOTJEILA aseYyg g" 
jusWI/4 
9 v A 0 CG p- 9 - 


e=bip#ibsy 4 
L=b'1#1611 





OOl- 


(QAne az 
$33139p) 
I 


0 


001 


00¢ 


87 


Single target, multiple harmonic cases in Section III.B (see 
Tables 10 and 13). Although the wavefront curvature is 
minimal at the lowest harmonic, another contribution to tie 
increase in range estimation error at the lower harmonics is 
the inefficiency introduced from the array design. 

From Tables 17 and 18, it appears that the [Moi 
algorithm is more sensitive to the inefficient array design 
at the lower harmonics. In the single harmonic test cases 
previously discussed (Section III.A), the maximum reliable 


range estimate determined by the LMS algorithm for the 0.0 


curve 
aB noise case waS approximately 1.0 * p. Whereas, for the 


Same conditions, the LMS algorithm could adequately 


flat 
determine a range estimate out to a range value of 3.0 * 9p. 
Given a scenario with an increased number of harmonics, this 
disparity in the maximum range estimate performance of the 
two algorithms is amplified at lower harmonics. Asa result, 
even though Target 1 exhibits Significant wavefront curvature, 
the LMS,,,,. algorithm does not efficiently analyze the Cage 
at the low harmonics (gq = 1,3). 

The wavefront curvature of Target 2 1s negligible at both 
of its harmonics, g = 2,5. The major cause for this is the 
depression angles. Target 2 is located at a depression 
angle 6, = 71°. This target, although positioned at the NF/7FE 
interface region, is located close to endfire. As a result, 


the phase variation across the Y axis is negligible (see 


Figure 3.9). The LMS,,... algorithm cannot correctly esti 
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the range of a target given this minimal wavefront curvature. 
At the higher harmonic (q=5), the LMS,,,, algorithm gives@a 
good range estimate for this target. 

Target 3, the target located in the FF region, radiates 
harmonics gq = 4,6. For the 0.0 dB case (refer to Table 17), 
IMS,,.ye Produces a good estimate at the highest harmonic q = 
6. Both algorithms perform poorly for this target at g = 4. 
From Figure 3.10 it can be seen that significant wavefront 
curvature is evident at these harmonics even with the target 
positioned in the far-field. As seen in our earlier results, 
this can be attributed to the near broadside orientation of 
the target relative to the array (0, = 2°). The superior 


ie) 


performance of LMS at this harmonic, can be attributed to 


curve 
Significant wavefront curvature coupled with the fact that 
the array design is optimized at this highest harmonic. The 
IMS,,,. algorithm adequately estimates the target rangewgae 
harmonic 6. In summary, this case is the most complex to 
analyze in that several factors must be considered simul- 
taneously at each harmonic to evaluate algorithm performance. 
At low harmonics, the performance of the LMS,,,, algorithm is 


Slightly superior to that of LMS for the given targets. 


curve 


The LMS,,.,. algorithm performed better than LMS). at uae 
highest harmonic for the target near broadside. Target 2 
generates minimal wavefront curvature and was most difficult 


to localize in range. All angular estimates were within five 


degrees for the 0.0 dB noise case and 10 degrees for -3.0 dB. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The purpose of this thesis was to develop a localization 
algorithm that could process spherical wave information 
incident upon a planar array and estimate the three spherical 
coordinates range r,, depression angle @,, and bearing angle 
ee This goal has been achieved by the signal processing 
model described in Section II.B and the localization algorithm 
outlined in Section II.C. Two different modified frequency- 
domain LMS adaptive algorithms were used to cancel the noise 
corruption of the output electrical signals. A comparison of 
the performance of each of these adaptive algorithms was 
detailed in the results section (Section III). The following 
comments highlight the significant findings from our results: 

- The spherical wave propagation model used in generating 
the ouput electrical signals at each element in a planar 
array has been validated. 

- The adaptive beamforming and non-linear least-squares 
estimation routines used in the signal processing section 
have been tested satisfactorily in a no-noise and noise 
environment. 

* In the absence of noise, or in high SNR cases, the 
localization models using both modified LMS adaptive 
algorithms accurately estimate range well into the far- 
field region. 

- Regardless of range estimation performance, both 
algorithms accurately estimated the bearing and depression 
angles in all cases tested. 

- In a noisy environment, for a signal target, the routine 
that assumed non-separable phase weights in the LMS 
algorithm performed best given a wavefront with 


Significant wavefront curvature. 


O:2 


*- In a noisy environment, for a single target, the routine 
that assumed separable phase weights in the LMS algorithm 
performed best given minimal wavefront curvature. Asa 
result, this routine could estimate the range to targets 
positioned at further distances from the receive array. 
Also, this routine was superior in estimating the range 
value at lower harmonics. Range estimates into the far- 
field were accurately determined by this algorithm in a 
noise environment. 

¢- Aside from the case where the target is located at an 
endfire position relative to the receive array, the full 
angular coverage capability of the algorithm was tested 
satisfactorily. 

- The ability to process several harmonics over the 
frequency spectrum of a single, broadband target was 
tested satisfactorily. 

* Given multiple targets positioned in both the near-field 
and far-field with unique spectral lines, both algorithms 
determined adequate estimates of the spherical coordinates 
for these targets in a noise environment. 

Given the design limitations presented in the introduction 
to Section III, the following trend was observed in the 
maximum range estimate capability of the separable weight 
routine LMS,,,,: In a noise environment characterized by a SNR 
value of 0.0 dB, the range estimate is valid to a distance of 
three times the range to the near-field/far-field boundary. 
For an 11 xX 11 element array designed to process a maximum 
frequency of 1000 Hz, the planar array dimensions would be 7 
meters x 7 meters. The maximum range estimate, as specified 
by the rule above, works out to be 180 meters (about 25 times 
the array length). If we use typical dimensions for a towed 


array, the maximum range value for this scenario would equal 


approximately 5 kilometers. 


= 


This maximum range rule for 0.0 dB SNR must be considered 
in light of the design limitations of our computer simulation 
(see Section III). Assuming that we have a computer system 
with more available resources (CPU's), we could increase two 
significant values: the number of time samples taken per 
element over the data record length, and the number of 
elements in the receive array. Our research has shown that 
by increasing the number of samples, our location estimates 
improve. With more elements, the array would be larger, 
thereby improving the detection of wavefront curvature. 
Additionally, with more elements, the performance of the 
Signal processing model is enhanced. One direct result of 
increasing these two parameters is an increase in the maximum 
estimation range determined by the algorithm. Our research 
has shown that the range estimate in a noise environment 
approaches the ideal no-noise result as we increase these 
parameters. Each test case conducted was validated to a range 
in excess of 100 kilometers for the no-noise condition. 

In the course of our analysis the following topics for 
further research surfaced: 

¢ Analyzing the performance of the localization algorithm 
in estimating target coordinates in an inhomogeneous 
medium. 

- Development of a surface fitting routine to replace the 
non-linear least-squares estimation algorithm which only 
performs smoothing along the X and Y axes. 

- Improving the performance of the modified LMS frequency- 


domain adaptive algorithms or replacing these algorithms 
with other beamforming algorithms. 
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